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ABSTRACT 

The  dynamic  evolution  of  tracked  waves  by  a  front-tracking  algorithm 
may  lead  on  either  numerical  or  physical  grounds  to  intersections  of  the 
waves.  The  correct  resolution  of  these  intersections  is  described  locally  by 
the  solution  of  Riemann  problems  and  requires  a  bifurcation  of  the  topology 
defined  by  the  tracked  waves. 

Here  we  describe  an  algorithm  that  is  appropriate  for  the  resolution  of 
scalar  tracked  waves,  such  as  material  discontinuities,  contact  discontinuities 
in  gas  dynamics,  or  constituent  concentration  waves  including  oil-water  banks 
in  oil  reservoirs.  Even  here  the  algorithm  is  not  fully  general,  and  the  resolu- 
tion of  the  intersections  of  an  arbitrary  set  of  curves  in  the  plane  for  the 
above  range  of  physical  problems  remains  unsolved.  However  with  the 
assumption  that  the  set  of  intersections  to  be  resolved  is  a  small  perturbation 
(resulting  for  example  from  a  small  time  step  in  an  evolution)  of  a  valid, 
non-intersecting  front,  the  algorithm  seems  to  be  general.  In  any  case  exam- 
ples will  be  presented  that  show  that  complicated  interfaces  can  be  generated 
automatically  from  simple  ones  through  successive  bifurcations. 
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1.    Introduction 

Front  tracking  is  a  computational  method  based  on  nonlinear  wave  analysis.  The  non- 
linear waves  and  their  interactions  are  conceptual  building  blocks  which  serve  to  organize  the 
computation.  This  method  is  characterized  by  high  resolution  solutions  for  problems  with 
significant  solution  discontinuities,  such  as  shock  waves  and  material  boundaries.  In  addition 
to  its  computational  achievements,  listed  for  example  in  [1,2,3,4,5],  the  development  of 
front  tracking  has  lead  to  advances  in  the  mathematical  theory  of  nonlinear  wave  analysis  and 
to  the  solution  of  a  long  standing  problem  for  pure  finite  differences:  the  20  year  old  Rach- 
ford  problem  of  viscous  fingering  for  immiscible  displacement  with  capillary  diffusion.  [6] 
The  present  article  is  devoted  to  the  solution  of  a  problem  which  had  been  viewed  as  a  major 
stumbling  block  for  front  tracking,  namely  the  intersections  of  fronts  with  other  fronts  and 
with  boundaries.  The  purpose  of  this  article  is  to  provide  a  solution  valid  for  the  case  of 
tracked  scalar  waves.  Partial  results  for  vector  wave  (shock)  interactions  have  been  reported 
elsewhere.  [1] 

We  distinguish  between  bifurcations  in  the  computed  approximate  solution  that  arise  on 
numerical  grounds  and  do  not  have  a  counterpart  in  the  exact  mathematical  solution  of  the 
underlying  equations,  and  those  that  do  correspond  to  bifurcations  in  the  exact  mathematical 
solutions  of  the  underlying  equations.  In  the  first  category  are  material  interfaces  and  contact 
discontinuities.  Within  the  limits  of  the  governing  equations  (neglecting  higher  order  effects 
such  as  diffusion  and  surface  tension),  such  interfaces  should  come  arbitrarily  close  and  form 
filaments,  without  ever  intersecting.  Such  filaments  are  numerically  difficult  to  compute  and 
usually  of  little  importance,  so  that  an  approximate  solution  method  that  allows  intersections 
to  occur  and  then  resolves  them  is  acceptable.  From  a  mathematical  point  of  view,  this 
means  simply  that  convergence  is  not  sought  in  the  Z.^  norm  but  in  L^  instead.  The 
occurrence  of  these  intersections  may  be  associated  in  part  with  the  order  of  accuracy  of  the 
time  integration  or  with  the  size  of  the  time  step  and  in  part  with  the  nature  of  the  unstable 
modes  excited.  For  example  the  method  of  [7,8]  which  uses  a  high  order  time  integration 
scheme  almost  never  develops  actual  crossing  of  the  interfaces.  The  front  tracking  code  used 
in  our  calculations  will  develop  crossings  for  comparable  problems  in  which  unstable  fingers 
dominate,  though  such  crossings  are  relatively  few  in  number.  However  the  unstable  jets 
with  pinch  nodes  excited  reported  on  here,  lead  to  interface  crossings  on  a  very  regular  basis. 

The  other  case,  in  which  the  bifurcation  is  a  property  of  the  exact  mathematical  solu- 
tion, can  occur  when  the  tracked  discontinuity  curve  represents  a  shock  wave.  We  distinguish 
between  vector  shock  waves  such  as  shocks  in  the  acoustic  mode  of  compressible  fluids,  and 
scalar  shocks  such  as  those  representing  concentration  or  saturation  fronts  in  oil  reservoirs, 
adsorption  problems,  or  chemically  reacting  fluids.    The  algorithm  discussed  in  this  paper  is 
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appropriate  only  for  tracked  scalar  shock  waves.    A  preliminary  discussion  of  the  bifurcation 
of  vector  shock  waves  is  given  in    [1]. 

Modifications  of  the  front  may  thus  be  necessary  for  both  physical  and  numerical  rea- 
sons. Algorithms  which  accomplish  such  tasks  have  been  referred  to  as  "untangling  algo- 
rithms" by  the  authors.  To  avoid  this  problem,  most  front  tracking  codes  use  some  version 
of  the  volume  of  fluid  method.  In  the  volume  of  fluid  codes  the  front  is  not  represented  by 
continuous  curves,  as  in  our  code,  but  is  reconstructed  from  a  marker  function,  which 
represents  the  amount  of  each  fluid  in  each  grid  block.  The  resulting  description  of  the  front 
is  of  a  lower  quality  than  it  is  for  front  tracking  as  defined  by  us,  see  [9].  For  example  the 
volume  of  fluid  method  may  cause  numerical  entrainment  and  fluid  mixing.  See  [10,11]  for 
a  further  discussion  of  interface  methods  and  see  [9,3,4,12,5]  for  additional  discussion  of 
front  tracking  as  used  in  this  paper. 

Over  several  years  we  have  experimented  with  preliminary  versions  of  untangle  algo- 
rithms which  have  lead  to  the  one  described  here.  A  related  untangling  problem  in  our  code 
occurs  in  the  generation  of  a  finite  element  grid,  used  to  solve  elliptic  equations.  An  initially 
rectangular  grid  is  locally  distorted  by  shifting  some  grid  points  along  grid  lines  to  adjacent 
front  curves  such  that  on  the  resulting  grid  the  front  lies  solely  along  the  boundaries  of  the 
finite  elements  [13].  In  a  complex  front,  extremely  close  curves  cause  an  ambiguity  for  the 
grid  point  shifting.  In  such  cases  both  the  curves  and  grid  are  locally  distorted  to  intersect  at 
a  common  position.  Here  the  curves  intersect,  although  they  do  not  cross,  and  these  curve 
crossings  have  to  be  resolved  before  the  elliptic  solver  is  invoked.  The  algorithm  to  accom- 
plish this  is  considerably  simpler  than  the  one  discussed  here,  and  has  been  described  in    [5]. 

2.   The  Algorithm 

To  establish  the  rules  governing  how  an  interface  is  to  be  untangled  it  is  necessary  to 
make  some  assumptions  about  the  tangled  interface.  Before  we  state  explicitly  the  restric- 
tions we  expect  to  hold  for  the  evolution  of  the  interface  let  us  establish  our  terminology. 
(See  [11]  for  a  more  complete  description  of  the  interface.)  We  restrict  our  discussion  here 
to  two  spatial  dimensions  though  the  restrictions  placed  on  the  front  evolution  are  equally 
valid  in  three  dimensions.  By  a  curve,  we  mean  a  piecewise  linear,  oriented,  curve  segment 
in  /?'.  The  linear  pieces  are  called  bonds.  The  start  and  end  points  of  a  curve  (which  need 
not  be  distinct)  are  called  nodes.  In  general  several  curves  will  start  or  end  at  a  common 
node.  The  incoming  curves  for  a  given  node  are  those  that  end  at  that  node,  while  the  outgo- 
ing curves  are  those  that  start  there.  An  interface  is  a  collection  of  curves  and  nodes.  A 
proper,  or  untangled  interface  is  one  whose  curves  intersect  (including  self-intersections)  only 
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at  their  start  or  end  points,  the  nodes.  Otherwise  the  interface  is  said  to  be  tangled.  These 
points  of  interior  curve  intersection  are  called  cross  points.  We  consider  two  interfaces,  an 
untangled  interface  l^  given  at  the  beginning  of  a  time  step  and  a  possibly  tangled  interface  /' 
defined  at  the  end  of  the  time  step.  The  interface  /'  is  connected  to  the  interface  /°  by  a 
linear  homotopy,  with  parameter  /,  0  <  r  <  1.  We  call  this  homotopy  the  time  step  homo- 
topy.  The  interfaces  1°  and  /'  separately  divide  the  plane  R^  into  regions  called  components. 
The  same  is  true  for  each  of  the  intermediate  interfaces  /'.  If  a  given  component  of  the  /' 
interface  does  not  connect  continuously  to  a  component  of  the  I^  interface  it  is  called  an 
unphysical  component.    Otherwise  it  is  physical. 

The  main  assumption  of  the  untangling  algorithm  described  here  is  that  the  curves  do 
not  cross  each  other  in  such  a  way  that  the  resulting  "tangled"  interface  /'  is  physically  mean- 
ingful. This  assumption  simplifies  the  identification  of  interface  tangling  having  occurred 
under  the  homotopy  by  ensuring  that  the  tangle  is  captured  in  the  form  of  crossed  curves  in 
the  /'  interface.    Our  explicit  assumptions  are: 

Al.  Each  point  on  a  curve  in  /'  has  crossed  some  other  curve  at  most  once  during  the 
last  time  step  in  which  /"  is  displaced  to  /'. 

A2.  Each  unphysical  component  has  on  its  boundary  at  least  one  cross  point  where  cer- 
tain curves  of  /'  intersect. 

A3.  At  the  cross  points  the  curves  really  cross  each  other,  rather  than  merely  touch. 

A4.  No  nodes  have  crossed  a  curve  during  the  time  step. 

Assumption  Al  is  really  a  limit  on  the  size  of  a  time  step.  A  step  that  causes  violation 
of  assumption  Al  is  too  big  and  must  be  reduced.  The  restrictions  implied  by  the  other 
assumptions  can  be  enforced  by  suitable  preprocessing  before  the  actual  untangling  algorithm 
is  started.  Assumption  A2  excludes  the  possibility  that  a  closed  loop  "inverts"  itself.  (See 
Fig.  1.)  Such  loops  can  be  identified  by  comparing  loop  orientation  at  t  =  0  and  t  =  1. 
Assumption  A3  is  included  for  completeness;  the  probabilities  of  a  point  lying  exactly  on 
another  curve  are  obviously  small.  Nodes  that  have  crossed  a  curve  must  also  be  identified 
in  advance  and  a  separate  algorithm  run  first  to  untangle  these  node  interactions.  Preprocess- 
ing of  this  sort  has  been  carried  out  to  some  degree  in  our  computations. 

The  untangling  procedure  then  consists  of  finding  the  cross  points  of  the  /'  interface, 
and  identifying  the  unphysical  areas  between  curves  that  have  crossed  each  other.  The  first 
step  in  the  algorithm  is  to  find  the  intersections  or  cross  points.  This  is  done  in  an  efficient 
way  described  in  [11].  Next  the  cross  points  are  replaced  by  nodes.  The  resulting  nodes  are 
called  non-homotopic  (NH)  nodes,  to  distinguish  them  from  the  (homotopic)  nodes  of  the 
interface  /',  which  are  connected  by  the  homotopy  to  the  nodes  of  the  interface  l^.    The  (old) 


curves  of  /'  that  cross  at  these  points  are  split,  and  the  new  curves  resulting  from  the  split  are 
attached  to  the  NH  nodes.  The  main  body  of  the  untangling  algorithm  then  consists  of  loop- 
ing over  these  NH  nodes,  and  determining  which  connected  curves  enclose  an  unphysical 
component.  In  doing  so  we  are  aided  by  several  observations  that  are  seen  to  follow  immedi- 
ately from  the  assumptions  Al  to  A4. 

01.  Every  NH  node  must  have  two  incoming  curves  and  two  outgoing  curves,  dividing 
the  area  surrounding  the  NH  node  into  four  angular  sectors.  An  incoming  and  outgoing 
curve  can,  however,  be  the  same  curve. 

02.  One,  and  only  one,  of  the  four  sectors  surrounding  an  NH  node  is  unphysical. 

03.  If  a  curve  starts  and  ends  at  the  same  NH  node  (i.e.  is  a  loop),  the  two  sectors  out- 
side the  curve  must  be  physical. 

04.  The  sector  between  two  curves  that  each  join  a  single  NH  node  n^  to  distinct  nodes 
n,  ^  n^,  must  be  physical  if  at  least  one  of  the  nodes  ^2  or  "3  is  a  homotopic  node. 

05.  An  unphysical  sector  must  be  bounded  by  either: 

(a)  a  single  curve  that  starts  and  ends  at  the  same  NH  node. 

(b)  two  curves  that  go  from  a  common  node  n^  to  another  common  node  rij,  one  of  which 
must  be  a  NH  node. 

We  will  not  prove  in  detail  any  of  these  observations,  although  they  are  perfectly  gen- 
eral conclusions  from  our  assumptions.  Observation  Ol  is  trivial,  since  a  cross  consists  of 
two  bonds  belonging  either  to  different  curves  or  the  same  curve  and  after  an  NH  node  has 
been  created  both  bonds  are  split  and  attached  to  the  node.  Observation  02  is  valid  from 
restrictions  Al  and  A3.  Observation  03  follows  from  observation  02,  since  for  a  curve  that 
starts  and  ends  at  the  same  NH  node,  the  component  exterior  to  this  loop  forms  two  of  the 
sectors  adjacent  to  the  NH  node.  However,  from  observation  02  one,  and  only  one,  can  be 
unphysical.    Observation  04  is  a  consequence  of  assumption  A4. 

The  implementation  of  the  untangling  algorithm  consists  of  two  parts.  The  first  part 
identifies  those  physical  components  that  can  be  unambiguously  determined  as  such  from  the 
topology  of  the  interface  /'.  For  those  NH  nodes  that  then  have  only  one  remaining  unidenti- 
fied sector,  this  must  be  the  unphysical  one  and  the  curves  bounding  it  are  therefore  deleted. 
The  second  part  of  the  algorithm  deals  with  nodes  that  have  two  or  more  adjacent  areas  that 
could  possibly  be  unphysical,  but  can  not  be  identified  correctly  without  reference  to  the  evo- 
lution. Rather  than  appeal  to  the  homotopy  defined  by  the  evolution  of  the  time  step,  we 
handle  these  nodes  in  a  somewhat  ad  hoc  manner,  using  the  experimental  observation  that 
the  unphysical  loops  are  usually  very  small.  We  shall  argue  though  that  the  way  we  treat 
these  nodes  is  reasonable. 


-6  - 

We  now  describe  the  algorithm  in  more  detail.  A  list  of  the  NH  nodes  inserted  at  the 
cross  points  is  generated.  The  list  also  contains  the  (four)  curves  connected  to  each  node, 
arranged  in  angular  order,  and  an  index  for  each  of  the  sectors  at  the  node  bounded  by  the 
curves.  A  first  loop  over  the  NH  node  list  applies  observations  03  and  04  to  identify  physi- 
cal areas,  and  sets  the  index  for  such  areas.  For  those  nodes  that  have  three  "physical"  sec- 
tors the  fourth  must  be  unphysical;  on  a  second  loop  the  curves  bounding  that  sector  are 
deleted  from  the  interface  /'  (and  from  the  NH  list),  and  the  node  removed  from  the  NH 
node  list.  Since  the  same  sector  can  be  adjacent  to  two  NH  nodes  (observation  05b),  it  is 
possible  to  now  encounter  NH  nodes  in  the  list  with  only  two  curves,  which  now  must 
represent  physical  curves.  We  therefore  check  for  such  nodes  in  the  loop,  and  delete  them 
from  the  NH  node  list.  This  second  loop  is  implemented  recursively,  finishing  when  no  NH 
nodes  with  three  identified  physical  sectors  can  be  found.  There  are  generally  only  a  few 
cross  points  at  each  time  step,  so  efficiency  in  this  part  of  the  algorithm  is  not  an  important 
consideration.  This  completes  the  first  part  of  the  algorithm.  The  remaining  entries  in  the 
NH  node  list  now  have  two  or  more  sectors  that  satisfy  observation  05. 

In  the  second  part  of  the  algorithm  we  find  the  area  of  each  sector  that  has  not  been 
marked  as  physical  in  the  first  part  of  the  code.  We  then  rely  on  the  observation  (based 
upon  numerical  experiment)  that  the  sector  with  the  smallest  area  is  usually  the  unphysical 
sector.  The  application  of  this  observation  is  tempered  and  supplemented  by  topological  res- 
trictions such  as  the  requirement  that  the  left  and  right  components  of  the  remaining  sectors 
be  consistent  with  those  of  near  by  curves.  It  is  obviously  possible  to  improve  this  part  of  the 
algorithm  somewhat,  in  particular  by  making  sure  that  any  sector  that  is  identified  as  the  one 
with  smallest  area  at  a  particular  NH  node  also  be  the  sector  of  smallest  area  at  any  other  NH 
node  to  which  it  may  belong.  To  achieve  greater  improvement  in  this  part  of  the  algorithm, 
it  is  necessary  to  refer  to  the  topology  of  the  interface  7°,  Even  in  doing  so,  there  may  still  be 
cases  that  can  only  be  resolved  by  "common  sense".  Fig.  2  demonstrates  one  such  case.  The 
curve  marked  A  can  be  thought  of  as  a  vortex  sheet,  and  the  curve  marked  B  as  enclosing  a 
blob  of  different  material.  The  dotted  line  marks  the  position  of  B  with  respect  to  A  at  the 
next  time  step.  It  is  obviously  impossible,  by  topological  arguments  referring  only  to  /°  and 
/',  to  distinguish  between  paths  (a)  or  (b)  and  hence  whether  area  x  or  y  represents  the  phy- 
sically impossible  area.  Only  by  knowledge  of  the  entire  homotopy  0  <  t  <  I  (rather  than 
just  its  two  endpoints)  is  this  possible.  This  however  presents  a  burdensome  numerical  task 
in  general;  our  belief  is  that  the  smallest  area  criterion  is  a  reasonable  choice,  since  if  there  is 
any  doubt  about  which  loop  should  be  deleted,  the  best  procedure  is  to  repeat  the  time  step 
with  a  smaller  time  increment. 


In  the  problems  for  which  this  untangling  algorithm  has  been  used  the  interface  curves 
do  not  carry  essential  physical  variables  that  are  not  also  representable  by  the  interior  degrees 
of  freedom.  Thus  for  any  given  discontinuity  in  the  solution  wave,  there  is  an  option  as  to 
whether  it  will  be  tracked  or  not.  In  the  present  code,  we  have  chosen  to  delete  from  the 
tangled  /'  interface  the  curves  surrounding  the  unphysical  sectors  that  result  from  the  crossing 
of  one  curve  in  l^  with  another.  For  material  boundaries  between  distinct  fluids,  this 
amounts  to  the  deletion  of  thin  filaments.  For  some  problems  it  might  be  advantageous  or 
even  necessary  to  retain  these  filaments  as  a  single  tracked  curve.  (See  Fig.  3)  For  example 
if  the  curves  represent  vortex  sheets  in  an  inviscid  fluid  and  if  the  fluid  in  the  interior  (away 
from  the  curves)  is  constrained  to  be  irrotational,  then  a  single  curve  is  needed  to  represent 
the  net  vorticity  of  the  two  curves  that  have  collided.  The  modification  of  the  algorithm  to 
handle  such  cases  is  obviously  a  minor  one. 

3.    Examples 

The  purpose  of  the  untangling  algorithm  presented  here  is  to  extend  the  application  of 
front  tracking  codes  to  the  solution  of  "messy,  practical"  problems.  We  demonstrate  this 
capability  by  an  example  involving  repeated  tracked  wave  interactions.  The  examples  shown 
are  taken  from  a  study  of  the  performance  of  various  well  arrangements  for  secondary  oil 
recovery.  Figure  4  shows  the  evolution  of  the  oil-water  front  for  a  line  drive  geometry 
where  the  middle  injection  well  has  been  displaced  from  the  line  set  by  the  other  four  injec- 
tion wells.  The  set  of  runs  from  which  this  example  is  taken  were  made  to  investigate  the 
effect  of  such  a  perturbation  in  a  homogeneous  reservoir  for  various  viscosity  ratios.  When 
the  frontal  mobility  ratio  is  sufficiently  high  the  offset  of  the  well  is  amplified,  and  the  mid- 
dle well  breaks  through  before  the  others.  In  the  run  presented  here,  the  frontal  mobility 
ratio  M  (for  the  forward  displacement  of  oil  by  water)  satisfies  the  stability  bound  .M  <  1. 
The  injected  fluid  (water)  is  therefore  less  mobile  than  the  surrounding  fluid  (oil).  Although 
the  expanding  water  fronts  move  laterally,  and  eventually  merge,  they  do  not  completely 
block  the  flow  of  the  more  mobile  oil  from  behind  the  injection  wells,  which  is  drawn  into 
narrow  channels  between  the  water  banks.  The  results  suggest  that  these  channels  are 
unstable,  with  pinch  off  instability  as  the  predominant  mode,  in  conjunction  with  a  weaker 
meandering  "garden  hose"  instability. 

The  four  figures  shown  in  Fig.  4  are  the  positions  of  the  water  bank  at  times  ranging 
from  the  time  shortly  after  the  first  merging  of  fronts  from  the  different  injection  sites  has 
taken  place  (a)  through  the  time  when  the  middle  well  has  broken  through  (d).  These  figures 
were  selected  such  that  one  can  trace  the  evolution  of  oil  "bubbles"  as  they  form  (b),  move 
toward  the  leading  stable  front  (c),  and  break  through  (d),  at  which  time  (in  this  particular 


example)  new  "bubbles"  are  again  forming. 

Instabilities  of  this  type  for  a  finger  of  viscous  fluid  moving  through  a  fluid  of  different 
viscosity  have  been  observed  previously  in  the  numerical  study  of  a  flow  of  two  fluids  in  a 
Hele-Shaw  cell  [7,8],  and  are  also  seen  experimentally  (for  recent  experiments  see  e.g. 
[14]).  Phenomena  with  a  related  appearance  occur  in  the  instability  and  breakup  of  jets  in 
compressible  gas  dynamics. 

There  are  two  noteworthy  differences  between  the  model  studied  here,  and  the  Hele- 
Shaw  model.  In  the  present  simulation,  no  physical  length  scale  comparable  with  the  role 
played  by  surface  tension  in  [7,8],  has  been  introduced  to  break  the  self  similarity  of  the 
governing  equations  at  small  length  scales.  Therefore  the  smallest  resolved  distance  (in  this 
case  the  grid  spacing)  plays  the  role  of  a  stabilizing  mechanism  at  small  length  scales. 

Fig.  5  shows  the  results  of  a  4  well  line  drive  study  with  the  inter-well  distances 
arranged  to  conform  with  that  of  Fig.  4.  The  frontal  mobility  ratio  of  this  series  of  runs. is 
slightly  different  from  Fig.  4  (0.288  compared  with  .586  for  Fig.  4).  The  permeability  func- 
tions and  connate  data  for  this  run  were  taken  from  a  particular  reservoir  field.  The  three 
plots  show  the  flow  computed  at  the  same  time  value  under  three  mesh  regimes  a)  5x20 
hyperbolic  8x30  elliptic  b)  10x40  hyperbolic  15x60  elliptic  c)  20x80  hyperbolic  30x120  elliptic. 
Under  grid  refinement  we  indeed  find  that  while  the  overall  evolution  is  the  same,  the  bub- 
bles become  smaller  and  the  time  of  their  first  appearance  in  the  computation  occurs  later 
than  in  coarser  grid  calculations  (This  later  conclusion  is  not  evident  from  Fig.  5  but  can  be 
seen  by  comparing  earlier  times  in  the  calculations  on  the  three  meshes.) 

For  realistic  reservoir  parameters  and  grid  spacings  the  capillary  pressure  length  scale, 
which  breaks  scale  invariance,  can  be  less  than  one  grid  spacing  and  so  grid  regularization  of 
this  small  length  scale  phenomena  is  appropriate.  The  grid  spacing  also  sets  a  length  scale 
for  numerically  induced  heterogeneity,  to  provide  a  seed  for  initialization  of  the  instability. 
Thus  while  we  believe  that  the  correct  zero-size-mesh  limit  for  the  equation  with  no  physical 
cut-off  of  short  distance  behavior  is  an  infinitely  thin  oil  channel  separating  the  water  banks, 
we  feel  the  numerical  grid  regularization  gives  a  more  correct  representation  of  a  flow  that 
would  result  with  a  physical  small  scale  cut-off  and  heterogeneity:  namely  the  appearance  of 
pinch  off  and  "garden  hose"  instabilities  in  this  channel. 

The  second  difference  is  the  rarefaction  region  behind  the  shock  in  the  model  for 
immiscible  displacement  calculated  here,  as  compared  with  only  a  contact  surface  in  the 
Hele-Shaw  model.  The  bubbles  must  traverse  this  rarefaction  region  as  they  move  towards 
the  leading  water  bank.  The  composition  of  the  bubble  will  thus  change  during  this  traversal. 
Both  the   pinch  off  (bubble  formation)   and  the  meandering  of  the  jet  cause  a  significant 
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amount  of  backflow  and  hysteresis  in  the  region  of  the  jet  between  the  main  fingers.  Back- 
flow  and  hysteresis  in  turn,  lead  to  a  rarefaction  fan  ahead  of  the  water  bank  in  the  channel. 
In  Fig.  6  the  saturation  contours  for  the  last  frame  in  Fig.  4  are  plotted.  The  contour  lines 
depict  0.1  increments  in  the  saturation.  Saturation  values  below  0.707  result  from  this  back- 
flow  and  hysteresis.  Thus  we  see  that  the  two  dimensional  jet  instability  provides  a  mixing 
mechanism  which  should  be  studied  in  more  detail. 

Our  last  figures  are  taken  from  a  study  of  more  typical  reservoir  well  placements  and 
contains  31  injection  and  production  wells.  The  fronts  and  well  positions,  early  in  the  calcu- 
lations are  shown  in  Fig.  7(a)  and  the  position  of  the  oil-water  front  when  the  last  production 
well  has  broken  through  is  shown  in  Fig.  7(b).  The  untangling  algorithm  handles  the  merg- 
ing of  the  banks  formed  from  water  injected  from  the  various  wells,  as  well  as  the  compli- 
cated case  where  many  banks  break  through  at  the  same  well. 

4.    Conclusion 

A  list  of  six  main  problem  areas  for  front  tracking  has  been  previously  mentioned  [4]: 

1.  Achieving  second-order  accuracy  at  a  shock  without  post-shock  oscillations. 

2.  Avoiding  stringent  limits  on  the  time  step  arising  in  the  bits  and  pieces  of  zones  which 
are  crossed  by  a  tracked  front. 

3.  Properly  treating  slip  along  a  front. 

4.  Treating 

a)  highly  distorted  fronts  and 

b)  changes  of  the  topology  of  regions  bounded  by  fronts  from  simply  connected  to 
multiply  connected  regions. 

5.  Treating  collisions  and  intersections  of  fronts  with  other  fronts  and  with  boundaries. 

6.  Treating  the  disappearance  of  weakening  fronts  and  the  appearance  of  new  fronts  at 
boundaries  or  at  collisions  of  other  fronts. 

The  status  of  front  tracking  on  these  six  difficulties  is  now  as  follows: 

1.  This  goal  requires  an  upgrade  in  the  riemann  solver  and  in  the  front  coupling  to  the 
interior  to  be  second  order  accurate. 

2.  This  problem  has  been  solved,  but  the  related  problem  of  optimal  front-interior  cou- 
pling deserves  further  study. 

3.  4a.  These  problems  are  fully  solved. 

4b,  5.  The  examples  given  in  this  article  demonstrate  the  capability  of  the  untangling 
algorithm  presented  here  to  handle  complex  interface  interactions. 

6.  Work  is  in  progress  on  these  problems.  What  is  required  is  the  implementation  of 
shock  capturing  methods  and  two  dimensional  Riemann  problem  solvers  [15]. 


-  10  - 


References 

1.  B.  Bukiet,  C.  L.  Gardner,  J.  Glimm,  J.  Grove,  J.  Jones,  O.  McBryan,  R.  Menikoff, 
and  D.  H.  Sharp,  "Applications  of  Front  Tracking  to  Combustion,  Surface  Instabilities 
and  Two-Dimensional  Riemann  Problems,"  ARO  Conference  Proceedings. 

2.  James  Glimm  and  D.  H.  Sharp,  "Numerical  Analysis  and  the  Scientific  Method,"  DOE 
Research  and  Development  Report  DOEIERI03077-270,  1986. 

3.  J.  Glimm,  B.  Lindquist,  O.  McBryan,  and  L.  Padmanabhan,  "A  Front  Tracking  Reser- 
voir Simulator  I:  The  Water  Coning  Problem,"  in  Frontiers  in  Applied  Mathematics, 
vol.  1,  SIAM,  Philadelphia,  1983. 

4.  I-L.  Chern,  J.  Glimm,  O.  McBryan,  B.  Plohr,  and  S.  Yaniv,  "Front  Tracking  for  Gas 
Dynamics,"  J.  Camp.  Phys.,  vol.  62,  pp.  83-110,  1986. 

5.  J.  Glimm,  W.  B.  Lindquist,  O.  McBryan,  and  G.  Tryggvason,  "Sharp  and  Diffuse 
Fronts  in  Oil  Reservoirs:  Front  Tracking  and  Capillarity,"  SIAM,  Proc.  Math,  and 
Comp  .  Methods  in  Seismic  Exploration      and  Reservoir  Modelling,  Houston,  Jan,  1985. 

6.  M.  J.  King,  W.  B.  Lindquist,  and  L.  Reyna,  "Stability  of  Two  Dimensional  Immiscible 
Flow  to  Viscous  Fingering,"  DOE  Research  and  Development  Report  DOEIERI03077-244, 
March,  1985. 

7.  G.  Tryggvason  and  H.  Aref,  "Numerical  Experiments  on  Hele-Shaw  Flow  with  a  Sharp 
Interface,"/.  Fluid  Mech.,  vol.  154,  pp.  1-30,  1983. 

8.  G.  Tryggvason  and  H.  Aref,  "Finger  Interaction  Mechanisms  in  Stratified  Hele-Shaw 
Flow,"  J.  Fluid  Mech.,  vol.  154,  pp.  287-301,  1985. 

9.  J.  Glimm,  E.  Isaacson,  D.  Marchesin,  and  O.  McBryan,  "Front  Tracking  for  Hyper- 
bolic Systems,"  Adv.  in  Appl.  Math.,  vol.  2,  pp.  91-119,  1981. 

10.  J.  M.  Hyman,  "Numerical  Methods  for  Tracking  Interfaces,"  Physica,  vol.  12D,  pp. 
396-407. 

11.  J.  Glimm  and  O.  McBryan,  "A  Computational  Model  for  Interfaces,"  Adv.  Appl. 
Math.,  vol.  6,  pp.  422-435,  1985. 

12.  J.  Glimm,  C.  Klingenberg,  O.  McBryan,  B.  Plohr,  D.  Sharp,  and  S.  Yaniv,  "Front 
Tracking  and  Two  Dimensional  Riemann  Problems,"  Adv.  in  Appl.  Math.,  vol.  6,  pp. 
259-290,  1985. 

13.  O.  McBryan,  "Elliptic  and  Hyperbolic  Interface  Refinement  in  Two  Phase  Flow,"  in 
Boundary  and  Interior  Layers  -  Computational  and  Asymptotic  Methods,  ed.  J.  Miller, 
Boole  Press,  Dublin,  1980. 


- 11  - 

14.  J.  V.  Maher,  "Development  of  Viscous  Fingering  Patterns,"  Phys.  Rev.  Lett.,  vol.  54, 
pp.  1498-1501,  1985. 

15.  W.  B.  Lindquist,  "Construction  of  Solutions  for  Two  Dimensional  Riemann  Problems," 
Adv.  Hyp.  PDE's.,  (To  Appear). 


-12 


Fig.  1.    An  example  of  the  inversion  of  a  loop  during  propagation  over  a  single  time 
step. 
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(a) 


Fig.  2.  A  case  where  two  curves  cross  each  other  in  which  it  is  impossible  to  determine 
the  area  which  represents  the  the  unphysical  region  from  only  the  topology  of  the  inter- 
face at  r°  or  r'. 
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(a) 


(b: 


(c: 


Fig.  3.    a)  A  typical  intersection  of  two  curves  which  define  the  sides  of  a  thin  material 
channel,    b)  One  method  of  resolving  the  resultant  intersection.    Here  the  channel  is 
pinched  off.    This  may  result  in  convergence  to  an  exact  solution  in  an  Z. ,  norm  rather 
than  L^.    c)  A  second  method  of  resolution  where  the  narrow  channel  is  represented  by 
a  single  curve.    Method  b)  is  currently  implemented  in  the  front  tracking  code  described 
in  the  text. 
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(a)  Time  step  22. 


(b)  Time  step  29. 


(c)  Time  step  36. 


(d)  Time  step  44. 


Fig.  4.    An  example  of  the  complex  interface  interactions  that  can  be  handled  by  our 
untangling  procedure:  a  line  drive  well  configuration.    The  calculation  is  for  immiscible 
displacement  with  quadratic  permeability  functions  and  equal  viscosities  of  oil  and 
water.    The  frontal  mobility  ratio  for  water  displacing  oil  is  0.586,  the  far  field  mobility 
ratio  is  1.0.    The  grid  for  the  elliptic  equations  is  32  by  32  mesh  blocks,  for  the  hyper- 
bolic 20  by  20. 
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(a) 


(b; 


(c) 


Fig.  5.    A  mesh  refinement  calculation  of  the  unstable  oil  channel  that  forms  between 
the  merging  water  banks.    Meshes  used  -  a)  5x20  hyperbolic,  8x30  elliptic  b)  10x40 
hyperbolic,  15x60  elliptic  c)  20x80  hyperbolic,  30x120  elliptic 
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Fig.  6.    Saturation  contours  in  0.1  increments  for  Fig.  4(d).    Selected  saturation  con- 
tours are  labeled. 
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(a)  Time  step  0. 


(b)  Time  step  40. 
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(c)  Time  step  80. 


(d)  Time  step  240. 


Fig.  7.    Plots  of  the  fronts  for  a  more  complicated  well  configuration  consisting  of  19 
injecting  wells  (crossed  squares)  and  12  producing  wells  (open  squares).    The  frontal 
mobility  ratio  for  water  displacing  oil  is  1.33.    The  elliptic  grid  is  64  by  64  mesh  blocks, 
while  the  hyperbolic  grid  is  48  by  48  mesh  blocks. 
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Fig.  8.    Saturation  contour  plors  for  Fig.  7(c).    The  region  in  the  rectangle  is  blown  up 
in  the  ne.xt  figure. 
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Fig.  9.    Blowup  of  a  subregion  in  Fig  8.    The  s  =   .2  contour  results  from  backflow. 
There  is  no  numerical  diffusion  across  the  tracked  front.    The  tracked  front  is  drawn  as 
a  heavier  line  than  the  saturation  contours. 
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